import numpy as np
import numpy.linalg as nla
import matplotlib.pyplot as plt
import nibabel as nib
import os
%matplotlib inline
path = '/home/snarles/predator'
os.listdir(path)
#os.listdir(path + '/8631_5_1_pfile/coil_images')
bvecs = np.loadtxt(path + '/chris1_bvec.csv', delimiter = ' ').T
b0_inds = list(np.where(np.array([nla.norm(v) for v in bvecs]) == 0)[0])
b0_inds = b0_inds[2:]
#plt.scatter(np.array(b0_inds),np.ones(10))
#wm = nib.load(path + '/8631_2_1_wm_mask.nii.gz').get_data()
wm = np.ones((120, 120, 69))
np.shape(wm)
c1_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil18_ec.nii.gz').get_data()
c1_0 = c1_0[:, :, :, b0_inds]
c2_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil21_ec.nii.gz').get_data()
c2_0 = c2_0[:, :, :, b0_inds]
c3_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil28_ec.nii.gz').get_data()
c3_0 = c3_0[:, :, :, b0_inds]
c4_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil31_ec.nii.gz').get_data()
c4_0 = c4_0[:, :, :, b0_inds]
c1 = c1_0
c2 = c2_0
c3 = c3_0
c4 = c4_0
for i in range(10):
c1[:, :, :, i] = c1_0[:, :, :, i] * wm
c2[:, :, :, i] = c2_0[:, :, :, i] * wm
c3[:, :, :, i] = c3_0[:, :, :, i] * wm
c4[:, :, :, i] = c4_0[:, :, :, i] * wm
c1_mu = np.mean(c1, axis = 3)
c2_mu = np.mean(c2, axis = 3)
c3_mu = np.mean(c3, axis = 3)
c4_mu = np.mean(c4, axis = 3)
c1_0mu = np.mean(c1_0, axis = 3)
c2_0mu = np.mean(c2_0, axis = 3)
c3_0mu = np.mean(c3_0, axis = 3)
c4_0mu = np.mean(c4_0, axis = 3)
(np.min(c1_mu[:]), np.max(c1_mu[:])), (np.min(c2_mu[:]), np.max(c2_mu[:])), (np.min(c3_mu[:]), np.max(c3_mu[:]))
vmin_ = -1.0
vmax_ = 15.0
figsize_ = (3, 3)
z_inds = [10, 20, 30, 40, 50]
#w_inds = range(10)
#z_inds = [12, 14, 16, 18, 20]
figsize_2 = (len(z_inds) * figsize_[0], figsize_[1])
f_c1, axarr = plt.subplots(1, len(z_inds), sharex=True, sharey=True, figsize = figsize_2)
for j in range(len(z_inds)):
z_ind = z_inds[j]
axarr[j].imshow(c1_mu[:, :, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
f_c2, axarr = plt.subplots(1, len(z_inds), sharex=True, sharey=True, figsize = figsize_2)
for j in range(len(z_inds)):
z_ind = z_inds[j]
axarr[j].imshow(c2_mu[:, :, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
f_c3, axarr = plt.subplots(1, len(z_inds), sharex=True, sharey=True, figsize = figsize_2)
for j in range(len(z_inds)):
z_ind = z_inds[j]
axarr[j].imshow(c3_mu[:, :, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
f_c4, axarr = plt.subplots(1, len(z_inds), sharex=True, sharey=True, figsize = figsize_2)
for j in range(len(z_inds)):
z_ind = z_inds[j]
axarr[j].imshow(c4_mu[:, :, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
w_inds = range(10)
vmin_2 = -1.0
vmax_2 = 1.0
figsize_2 = (len(w_inds) * figsize_[1]/2, 4 * figsize_[0])
f_sv, axarr = plt.subplots(4, len(w_inds)+1, sharex=True, sharey=True, figsize = figsize_2)
z_ind = 20
axarr[0, 0].imshow(c1_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
axarr[1, 0].imshow(c2_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
axarr[2, 0].imshow(c3_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
axarr[3, 0].imshow(c4_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
for i in range(len(w_inds)):
w_ind = w_inds[i]
axarr[0, i+1].imshow(c1[:, 20:80, z_ind, w_ind] - c1_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
axarr[1, i+1].imshow(c2[:, 20:80, z_ind, w_ind] - c2_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
axarr[2, i+1].imshow(c3[:, 20:80, z_ind, w_ind] - c3_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
axarr[3, i+1].imshow(c4[:, 20:80, z_ind, w_ind] - c4_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
z_inds = range(18, 26)
w_inds = range(10)
figsize_2 = (len(z_inds) * figsize_[0]/2, (len(w_inds)+1) * figsize_[1])
f_c1, axarr = plt.subplots(len(w_inds)+1, len(z_inds), sharex=True, sharey=True, figsize = figsize_2)
for j in range(len(z_inds)):
z_ind = z_inds[j]
axarr[0, j].imshow(c3_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
for i in range(len(w_inds)):
for j in range(len(z_inds)):
w_ind = w_inds[i]
z_ind = z_inds[j]
axarr[i+1, j].imshow(c3[:, 20:80, z_ind, w_ind] - c3_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
z_inds = range(18, 26)
figsize_2 = (len(z_inds) * figsize_[0]/2, 4 * figsize_[1])
f_c1, axarr = plt.subplots(4, len(z_inds), sharex=True, sharey=True, figsize = figsize_2)
vmin_3 = -2.0
vmax_3 = 2.0
for j in range(len(z_inds)):
z_ind = z_inds[j]
axarr[0, j].imshow((c3_0mu[10:100, 20:80, z_ind]) * wm[10:100, 20:80, z_ind],
cmap='bone', vmin = vmin_, vmax = vmax_)
axarr[1, j].imshow((c3_0mu[10:100, 20:80, z_ind] - c3_0mu[9:99, 20:80, z_ind]) * wm[10:100, 20:80, z_ind],
cmap='bone', vmin = vmin_3, vmax = vmax_3)
axarr[2, j].imshow((c3_0mu[10:100, 20:80, z_ind] - c3_0mu[10:100, 19:79, z_ind]) * wm[10:100, 20:80, z_ind],
cmap='bone', vmin = vmin_3, vmax = vmax_3)
axarr[3, j].imshow((c3_0mu[10:100, 20:80, z_ind] - c3_0mu[10:100, 20:80, z_ind-1]) * wm[10:100, 20:80, z_ind],
cmap='bone', vmin = vmin_3, vmax = vmax_3)
b0_breaks = np.array([2] + b0_inds)
break_starts = b0_breaks[:-1]+1
break_ends = b0_breaks[1:]-1
break_lens = break_ends - break_starts
len(break_lens)
def diffusion_aggregate(data_0):
data_d = np.zeros((120, 120, 69, 10))
for i in range(10):
data_d[:, :, :, i] = np.mean(data_0[:, :, :, break_starts[i]:break_ends[i]], axis = 3) * wm
return data_d
def demean(data_0):
data_0 = np.array(data_0)
wl = np.shape(data_0)[3]
data_mu = np.mean(data_0, axis = 3)
for i in range(wl):
data_0[:, :, :, i] = data_0[:, :, :, i] - data_mu
return data_0
def da_dac(data_0):
data_da = diffusion_aggregate(data_0)
data_dac = demean(data_da)
return data_da, data_dac
c1_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil18_ec.nii.gz').get_data()
c1_da, c1_dac = da_dac(c1_0)
del c1_0
c2_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil21_ec.nii.gz').get_data()
c2_da, c2_dac = da_dac(c2_0)
del c2_0
c3_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil28_ec.nii.gz').get_data()
c3_da, c3_dac = da_dac(c3_0)
del c3_0
c4_0 = nib.load(path + '/8631_5_1_pfile/coil_images/8631_5_coil31_ec.nii.gz').get_data()
c4_da, c4_dac = da_dac(c4_0)
del c4_0
w_inds = range(10)
vmin_2 = -1.0
vmax_2 = 1.0
figsize_2 = (len(w_inds) * figsize_[1]/2, 4 * figsize_[0])
f_sv, axarr = plt.subplots(4, len(w_inds)+1, sharex=True, sharey=True, figsize = figsize_2)
z_ind = 20
axarr[0, 0].imshow(c1_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
axarr[1, 0].imshow(c2_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
axarr[2, 0].imshow(c3_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
axarr[3, 0].imshow(c4_mu[:, 20:80, z_ind], cmap='bone', vmin = vmin_, vmax = vmax_)
for i in range(len(w_inds)):
w_ind = w_inds[i]
axarr[0, i+1].imshow(c1_dac[:, 20:80, z_ind, w_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
axarr[1, i+1].imshow(c2_dac[:, 20:80, z_ind, w_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
axarr[2, i+1].imshow(c3_dac[:, 20:80, z_ind, w_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)
axarr[3, i+1].imshow(c4_dac[:, 20:80, z_ind, w_ind], cmap='bone', vmin = vmin_2, vmax = vmax_2)